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ABSTRACT 

We have measured the infrared transit of the extrasolar planet HD 209458 b 
using the Spitzer Space Telescope. We observed two primary eclipse events (one 
partial and one complete transit) using the 24 pm array of the Multiband Imaging 
Photometer for Spitzer (MIPS). We analyzed a total of 2392 individual images 
(10-second integrations) of the planetary system, recorded before, during, and 
after transit. We perform optimal photometry on the images and use the local 
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zodiacal light as a short-term flux reference. At this long wavelength, the tran- 
sit curve has a simple box-like shape, allowing robust solutions for the stellar 
and planetary radii independent of stellar limb darkening, which is negligible at 
24 jum. We derive a stellar radius of R* = 1.06 ± 0.07 R©, a planetary radius 
of R p = 1.26 ± 0.08 Rj, and a stellar mass of 1.17 M©. Within the errors, our 
results agree with the measurements at visible wavelengths. The 24 /t: m radius 
of the planet therefore does not differ significantly compared to the visible re- 
sult. We point out the potential for deriving extrasolar transiting planet radii to 
high accuracy using transit photometry at slightly shorter IR wavelengths where 
greater photometric precision is possible. 

Subject headings: extrasolar planets, stars: individual (HD 209458) 


1. Introduction 

The transit of an extrasolar planet across its star allows us to measure the radius of 
the planet (Charbonneau et al. 2000; Henry et al. 2000; Brown et al. 2001). Of the ten 
known transiting planets, HD 209458 b, with a radius (at visible wavelengths) of 1.320 ± 
0.025 Rj (Knutson et al. 2006), is inflated compared to the other known transiting planets 
and thus has a lower bulk density. One explanation for the anomalous radius was inflation 
by the dissipation of tidal stress within the planet (Bodenheimer et al. 2001). However, the 
timing of the secondary eclipse as observed by the Spitzer Space Telescope (Deming et al. 
2005b), as well as improved radial velocity observations (Laughlin et al. 2005), have ruled 
out a non-zero orbital eccentricity of the magnitude (~0.03) needed by the tidal dissipation 
theory. Showman & Guillot (2002) suggest that kinetic energy produced by atmospheric 
circulation and deposited in the planet’s interior could account for the missing energy source 
and increase the planetary radius. Another popular proposed explanation is the possibility 
of obliquity tides (Winn & Holman 2005), in which a non-zero obliquity (made possible by 
a spin-orbit resonance) could drive the tidal dissipation and provide the necessary energy to 
give the planet an inflated radius. 

Our Spitzer program to measure the 24 pm flux of HD 209458 b includes observations 
during transit (i.e., primary eclipse), revealing the infrared (IR) radius of the planet, which 


1 Thi$ work is based on observations made with the Spitzer Space Telescope, which is operated by the 
Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this 
work was provided by NASA. 
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is reported in this paper. Section 2 further elaborates on why an IR radius measurement 
is of interest. Section 3 describes the observations; section 4 explains the photometric data 
analysis and radius fit. Section 5 concludes with results and discussion. 


2. Motivation for Infrared Radius Measurements 

The reduced stellar flux at mid-IR wavelengths implies that transit photometry in this 
region is unable to achieve the high photometric precision obtained at visible wavelengths 
(Brown et al. 2001). However, stellar limb darkening weakens with increasing wavelength 
due to the increasing H~ free-free opacity (Vernazza et al. 1976). Thus, the fitting of transit 
curves to mid-IR data is simple and robust, and gives results independent of limb-darkening 
parameterizations. A mid-IR radius measurement is also of intrinsic interest for understand- 
ing the planet. Observations of the planet during transit suggest that clouds and scattering 
layers could potentially extend to great heights. The low sodium abundance reported by 
Charbonneau et al. (2002) and the upper limit on CO reported by Deming et al. (2005a) 
both support this physical picture of the planet’s atmosphere. Furthermore, the escaping 
atmosphere observed by Vidal-Madjar et al. (2003) suggests that other processes may be at 
work in planet’s atmosphere. The scale height of the atmosphere is ~450 km, and if the 
clouds extend to even a few scale heights, this would represent ~0.01 Rj. While our obser- 
vations are not able to reach that level of precision in the planetary radius, it is nonetheless 
valuable to search for unexpectedly large variations in the planetary radius as a function 
of wavelength. A major goal of our observational program was to test this scenario by 
measuring the 24 fii n radius of HD 209458 b. 


3. Observations 

We observed two transit events using the Spitzer Space Telescope (Werner et al. 2004): 
a half eclipse event (ingress only) on 2004 December 5 and a full eclipse event on 2005 
June 27. We used the MIPS (Rieke et al. 2004) 24 yum array, which is a Si: As detector 
with 128x128 pixels and an image scale of 2.55 arcsec pixel -1 . We obtained two series of 
10-second exposures using the standard MIPS raster pattern, which places the star at 14 
different positions on the array. This produced 864 images during the half eclipse event and 
1728 images during the full eclipse event. 
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4. Analysis 
4.1. Photometry 

We first reject obviously bad images, including those with strong ‘jailbar’ features 2 , as 
well as the initial image of each cycle due to a prominent ‘first frame’ effect. This leaves 780 
images for the half eclipse event, and 1612 images for the full eclipse event. We analyze each 
eclipse event separately. For each of the 14 raster positions, we perform the following steps: 

1. Median filter the images to remove energetic particle events and hot pixels. The median 
image is computed from all the images (typically 60) at a given raster position, and the 
difference image is constructed by subtracting the median image from a given image. 
We then employ a routine called SIGMAJ7ILTER 3 , which computes the mean and 
standard deviation of the pixels within a box of specified width (excluding the center 
pixel). If the center pixel is deviant by more than a specified number of standard 
deviations from its neighbors, it is replaced by the mean of the remaining pixels in the 
box. We use a box width of 20 pixels and a sigma limit of 10, and iterate until no 
more pixels are changed. This cleans most of the hot pixels evident in a given image. 
Rejected from further consideration are any images in which a pixel within the defined 
aperture containing the star is changed. 

2. Calculate and subtract the total zodiacal background from all pixels in a given image. 
The background level for each image is determined by constructing a histogram of 
all pixels in an image (with a bin size of 0.01 MJy/sr) and fitting a Gaussian to the 
histogram. The center of the Gaussian fit is then the ‘average’ background level for 
that image, and this constant level is subtracted from the each pixel in the image to 
create the background-subtracted image. 

3. Find the center of the star to a precision of 0.01 pixel by dithering the theoretical 
point spread function (PSF) over the individual images and finding the best fit. The 
theoretical PSF was obtained 4 for a 5000 K blackbody on the center of the 24 pm 
array. These files have been modeled using Tiny Tim (Krist 1993). We first resample 
the theoretical PSF to a scale 100 times finer for comparison to the real data. The 
resampled PSF is dithered in both dimensions, and the ‘best-fit’ PSF to the data (using 
linear least squares) determines the center of the star. 

2 See the MIPS Data Handbook, available at http://ssc.spitzer.caltech.edu/mips/dli/ 

3 See the IDL Astronomy Library at http://idlastro.gsfc.nasa.gov/contents.html 

4 See http://ssc.spitzer.caltech.edu/mips/psffits/ and the file mips_24.5000c .fits. 
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4. Use the best-fit PSF to weight the pixels near the star before adding them. This is 
applied to the background-subtracted images and produces optimal photometry (anal- 
ogous to Horne 1986). The errors are derived by propagating the MIPS errors through 
the optimal error formula (Horne 1986, see his Table 1, item 8). The errors are domi- 
nated by statistical fluctuations in the zodiacal background. The optimal photometry 
typically provides a S/N improvement of 50% or more over the standard aperture 
photometry approach (where the pixels are summed with no weighting function). 

5. Normalize the optimal photometry to the total background level in the frame. The total 
background level is simply the average of all the points (except those in a 3x3 box 
surrounding the star) in a given image. This gives the stellar intensity relative to the 
zodi and thus removes any remaining instrument response variations. 

Having completed these steps for each raster position, we recombine the data to obtain the 
entire time series for each eclipse event. Time is converted to orbital phase using the most 
recent and most accurate orbital period and ephemeris (Knutson et al. 2006), and we account 
for the light travel time between the Sun and Spitzer. The zodiacal background changes 
linearly by ~1.5% over the six-hour duration over which the full eclipse event is observed, 
and we remove this effect from both eclipses separately. We estimated the magnitude of this 
change using the Spitzer Observations Planning Tool (SPOT), which employs the zodiacal 
dust model from Kelsall et al. (1998). Finally, the time series are normalized to an average 
of unity for the out-of-transit points. 

The calibrated, unbinned photometry is given in Table 1 and is shown in Figure 1 (upper 
panel). The upper panel shows the aggregate data for both events combined (2392 points), 
and it clearly reveals the eclipse. The lower panel shows the average in bins of phase width 
0.001; the box-like shape of the light curve due to the lack of stellar limb darkening is quite 
evident. To perform the optimal photometry, we used two independent codes, derived from 
the same basic algorithm but constructed by individual researchers. We obtained virtually 
identical per-point results with both. 


4.2. Light Curve Fitting 

We construct a family of simple, approximate light curves by connecting intensities at 
the contact times with straight line segments. Four observable parameters uniquely describe 
a light curve in the absence of stellar limb darkening: 1) the duration of full eclipse (i.e., 
the time between second and third contact, ip), 2) the total duration of the eclipse (i.e., 
the time between first and fourth contact, £t)> 3) the eclipse depth, and 4) the observed 
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time of center eclipse. We derive best-fit values for these four observables by minimizing 
the reduced chi-squared (xl) of the fits to each simple transit curve generated. We refer 
to this ‘trial-and-error’ technique as ‘Method 1.’ We have chosen a large enough range of 
parameters and small enough grid spacing to avoid finding only a local minimum. Figure 2 
shows xl & s a function of tp and tp for zero phase offset (the best-fit value); the upper panel 
shows the contour plot for the xt, and the lower panel casts the result in terms of confidence 
intervals. (The contours and confidence intervals are correctly calculated by projecting the 
xl into the plane of interest; see Press et al. (1992, Section 15.6 and Figure 15.6.4).) The 
best- fit observables (tp, tp, eclipse depth, and time offset) are listed in Table 2 (column 
marked ‘Method 1’). 

In order to verify the results from Method 1 and to ensure that we did not hit a local 
minimum, we employed a second method of finding the best-fit observables to the data, 
the MPFIT package 5 , which performs a Levenberg-Marquardt least-squares fit. Combined 
with our own function that computes a theoretical eclipse (as described above) , this method 
independently calculates the four observables, and these results are also shown in Table 2 
(column marked ‘Method 2’). We note that the two results agree closely, and we adopt the 
results from MPFIT (Method 2) as our formal results. Although the resulting minimum xl 
is slightly larger, we nonetheless adopt Method 2 because it is computationally more efficient 
for calculating the errors on the individual parameters, as discussed below. 

Using the analytic formulation of Seager <k Mallen-Ornelas (2003) , we derive the physical 
parameters of the system from the observables (tp, tp, and eclipse depth, obtained from 
the best-fit simple light curve, and the known orbital period (Knutson et al. 2006)). The 
impact parameter b, the ratio of the semi-major axis to the stellar radius, a/R x , the orbital 
inclination i, and the stellar density p* are derived immediately from these observables 
(Seager &: Mallen-Ornelas 2003, see their Equations 7, IT— T9) 6 . 

Determining the stellar and planetary radii from these parameters requires an assump- 
tion of the stellar mass (Brown et al. 2001), We assumed stellar masses between 0.9 and 
1.30 M 0 (shown in Figure 3), covering a region surrounding the reported stellar mass of 
1.146 M 0 (Brown et al. 2001). Interestingly, by using the analytic formulation, we note that 
the orbital inclination is determined by the transit times and the ratio of the stellar radius 
to the orbital semi-major axis; this means that our assumption of the stellar mass, while 
determining the stellar radius, does not affect the orbital inclination. 


5 http : //cow . physics . wise . edu/~craigm/idl/idl . html 

6 Note that their Equation 19, derived from Equation 9, is missing a factor of |tt. 
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We use a bootstrap Monte Carlo method (Press et al. 1992, see Section 15.6) to deter- 
mine realistic errors in the observables and physical parameters. We randomly select N=2392 
data points with replacement (meaning some points are duplicated) to create a ‘synthetic’ 
data set. We create 10000 such synthetic data sets, perform the same fitting procedure 
described above (again, Method 2, using MPFIT) to each one, and derive the physical pa- 
rameters from the best-fit observables. In this way, we derive a set of physical parameters 
(R*, R p , and i) for each of the 10000 synthetic data sets. Using the routine HISTOGAUSS 
(from the IDL Astronomy Library 7 ), we fit a Gaussian to each one of the arrays of observ- 
ables and physical parameters. All parameters are normally distributed and symmetric, so 
that the width of each best-fit Gaussian represents the 1-a error in the associated parameter, 
and these are the uncertainties presented in Tables 2 and 3. 

The bootstrap method was also performed on each of the eclipse events separately to 
determine the observed time of center eclipse, as shown in Table 4. Both are consistent with 
zero offset in time from the predicted value. For the half eclipse (event 1), we set the full 
eclipse time (£p) to the value derived from fitting the aggregate data and hold it fixed, while 
minimizing the other three observables. As expected, the uncertainty in the time of center 
for the half eclipse is much larger than that of the full eclipse. 

Next, we check the radii and orbital inclination by removing the approximation that the 
light curve is comprised of straight line segments. We have developed a routine to compute 
light curves numerically (Richardson et al. 2006), and we include the small effect due to 
predicted limb darkening at 24 pm, derived from a Kurucz 8 model atmosphere for stellar 
parameters T e = 6000 K, logy = 4.5, and [Fe/H]= 0.0. We validated the code by verifying 
that we can reproduce the fits to the very precise HST optical data from Brown et al. (2001). 
We adopted the derived parameters from the best-fit simple curve and calculated an exact 
theoretical light curve. The result is plotted as a dashed line in the lower panel of Figure 1, 
but it is hard to see since it is nearly identical to the simple curve. The xl °f the fits for each 
of the two curves to the data are nearly identical: 1.0060 for the best-fit simple curve from 
MPFIT, compared to 1.0061 for the theoretical light curve. We therefore conclude that the 
limb darkening is negligible at 24 /tm (as expected) and that the simple light curve composed 
of straight line segments is an accurate method of deriving the physical parameters. 

Finally, we checked our results by incorporating information from the transit at visible 
wavelengths (Brown et al. 2001). There the transit depth is 0.0164, compared to 0.0149 ± 


7 http : //idlastro . gsf c . nasa, gov/contents . html 

8 Available from http:/ /kurucz.harvard.edu. We linearly interpolate the Kurucz parameters at 20 and 
40 pm to estimate the values at 24 pm. 
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0.0003 at 24 /on, a significant difference. The ratio of visible to IE transit depth can be 
used to determine the impact parameter, the minimum projected radius where at which 
the planet crosses the star, and thereby the orbital inclination. That is, we are deriving 
the degree of limb darkening at the given projected stellar radius of closest approach and 
using the Kurucz model to determine the location of the chord. We used the limb darkening 
tabulated by a Kurucz model atmosphere (same parameters used above), and the ‘small 
planet’ approximation from Mandel & Agol (2002). We calculated numerically the ratio of 
visible to IR transit depth as a function of impact parameter. Comparing the observed ratio 
(1.100±0.03) to this relation gives an impact parameter of 0.58±0.07, and an orbit inclination 
of 86.6° ±0.6°. Within the errors, this agrees with the results at visible wavelengths (Brown 
et al. 2001; Wittenmyer et al. 2005) and is consistent with the i = 87.97° ± 0.85° value we 
derive internally from our IR data. This calculation serves as an independent check of our 
results and a direct comparison to the visible results. 


5. Results and Discussion 

We have computed the stellar density directly from the observable quantities from the 
best-fit simple curve (Seager & Mallen-Ornelas 2003). Assuming a stellar mass allows us to 
calculate the stellar radius. This empirical mass-radius relation is shown in Figure 3, where 
we have derived the radii for stellar masses from 0.9 to 1.3 M 0 . We break the degeneracy 
by intersecting the stellar radius curve with the mass-radius relation from Cody & Sasselov 
(2002), which is derived by fitting stellar models to a constant luminosity. This is shown 
as the dashed line in Figure 3, following Wittenmyer et al. (2005, their Figure 5). On this 
basis, we derive the stellar mass to be M = 1.171 M 0 , with R* = 1.06 ± 0.07 R 0 and R p = 
1.26 ± 0.08 Rj. Our result for the planetary radius agrees with the updated visible radius 
of R p = 1.320 ± 0.025 Rj (Knutson et al. 2006). 

We are encouraged by the fact that our radius error is only four times larger than that 
obtained by Knutson et al. (2006), in spite of the fact that our infrared photometric precision 
is an order of magnitude poorer than the HST visible photometry. We attribute this to the 
character of the IR transit curve, where the lack of limb darkening produces a simple transit 
shape, from which radius information can be extracted with maximum efficiency. We point 
out that photometry at other accessible Spitzer wavelengths such as 8 and 16 fim would 
provide much higher photometric precision for bright transiting systems, because the stellar 
flux will be much higher, and the zodiacal background will not be a limiting factor. Limb 
darkening remains sufficiently weak at these shorter IR wavelengths to maintain a relatively 
simple transit light curve shape. Considering also that Spitzer’s heliocentric orbit allows 
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uninterrupted observations of complete transits, we suggest that IR transit photometry from 
Spitzer may be the optimal method for precise radius determination in bright transiting 
planet systems. 


This work is based on observations made with the Spitzer Space Telescope, which is 
operated by the Jet Propulsion Laboratory, California Institute of Technology under a con- 
tract with NASA. Support for this work was provided by NASA. LJR acknowledges support 
as a NASA Postdoctoral Fellow (formerly NRC Research Associate) at NASA Goddard. We 
thank the referee for insightful comments and suggestions that significantly improved the 
manuscript. 

Facilities: Spitzer(MIPS) 
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Table 1. Calibrated, unbinned photometry. 


HJD 

Phase 

Relative Intensity 

Error 

2453344.6466211 

-0.03728 

1.009530 

0.008265 

2453344.6467485 

-0.03724 

1.001394 

0.008643 

2453344.6468759 

-0.03721 

1.010476 

0.008172 

2453344.6470034 

-0.03717 

1.012307 

0.008398 

2453344.6471308 

-0.03713 

0.992136 

0.008321 


[The complete version of this table is in the electronic edition of the Journal. The printed 

edition contains only a sample.] 
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Table 2. Observables derived from xl minimization of simple eclipse curves 


Parameter 

Method 1 

Method 2 

Error 

Depth 

0.01496 

0.01493 

0.00029 

t T (hr) 

2.978 

2.979 

0.051 

t F (hr) 

2.254 

2.253 

0.058 

Time Offset (hr) 

0.000 

0.001 

0.013 

xl 

1.00514 

1.00598 

- 
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Table 3. Derived physical parameters for the two minimization techniques. 


Parameter 

Method 1 

Method 2 

Error 

Assumed Stellar Mass (M 0 ) 

1.173 

1.171 


Stellar Radius (R 0 ) 

1.063 

1.064 

0.069 

Planetary Radius (Rj) 

1.265 

1.265 

0.085 

Orbital Inclination (deg) 

88.00 

87.97 

0.85 
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Table 4. Observed time of center eclipse. 



Time (HJD) 

Error 

Event 1 (Half) 

2453344.768245 

0.002608 

Event 2 (Full) 

2453549.201422 

0.000617 
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Fig. 1.— Upper Panel: All 2392 measurements versus heliocentric phase. Lower Panel: Data 
averaged in phase (bin size = 0.001 in phase); also shown are the best fit straight-line curve 
(solid line) and the exact theoretical light curve (dashed line, difficult to see), calculated using 
the best-fit physical parameters derived from the straight-line curve. Heliocentric phase was 
computed using the orbital period and ephemeris from Knutson et al. (2006). 
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Fig. 2. — Upper Panel: Contour plot of the reduced chi-squared (xt) fit from Method 1, 
showing the total eclipse time vs. the full eclipse time at zero phase offset (best-fit value). 
Lower Panel: Same result, but converted to confidence interval in standard deviations. In 
both panels, the minimum xt (1.0051) is marked by an X. 
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Fig. 3. — The stellar density derived from the best-fit parameters can be transformed to 
stellar radius by assumption of a range of stellar masses. The transit data therefore yield 
an empirical mass-radius relation. Upper Panel: Planetary radius as a function of assumed 
stellar mass. Solid line represents result from MPFIT (Method 2); dotted line represents 
the ‘trial-and-error’ minimization technique (Method 1). Lower Panel: Stellar radius as a 
function of assumed stellar mass. Methods 1 and 2 indicated as in upper panel. Dashed line 
represents mass-radius relation from Cody &; Sasselov (2002). Intersection of this relation 
with the empirical curve to the assumed stellar masses allows a determination of the stellar 
mass (1.171 M©), marked by the vertical dotted line. This reveals the best fit stellar and 
planetary radii, R* = 1.06 R© and R p = 1.26 Rj. 




